%% PDE SOLVER FOR THE KNWON-TYPE CONTRACT
function [V, C, Beta, K, Eta, f, ExRet, mu_v] = type_system_solver(params, grid_params, num_params)
%% PARAMTERS

%parameters
temp = num2cell(params);
[rho, r, ra, alpha, sigma, lambda, M, eta_min, eta_max] = temp{:};

[n_grid, Y_mat, PHI_mat, dy, dphi] = grid_params{:};

[max_iter, D, thr] = num_params{:};

PHI = PHI_mat(1,:)';

%% Initialize matrices
V_out = (((ra-r)*(1-rho)/rho)^(rho/(1-rho))); %Value for shutting down the firm
V00 = V_out*ones(n_grid, 1);

C0 = 0.01*ones(n_grid,1);
Beta0 = 0.01*ones(n_grid, 1);
Eta0 = eta_max*ones(n_grid, 1);

%% ITERATIONS

iter = 1;
reiterate = 1;
while reiterate
            
    tic        

    %% Find optimal policy.
    [C, Beta, Eta, K, f, ExRet, mu_v] = type_optimal_policy(V00, C0, Beta0, Eta0, params, grid_params);
    
    %% Find value function
    V = type_hjb_solver(V00, Beta, f, mu_v, params, grid_params, num_params);

    %% Check change
    vmin = min(V);
    vmax = max(V);
    
    [posmin] = find(V == vmin,1); 
    [posmax] = find(V == vmax,1); 
    
    disp(['Min of ', num2str(vmin), ' at ', num2str(posmin)])
    disp(['Max of ', num2str(vmax), ' at ', num2str(posmax)])
    
    disp(['Step size of ', num2str(D)])

    changemat = abs(V./V00 -1);
    test = max(changemat);
    [poschange] = find(changemat == test,1);
    disp(['Change of ', num2str(test,'%.4e'), ' at ', num2str(poschange)])

    if  test < thr
        reiterate = 0;
        disp(['Converged at iteration ',num2str(iter)])
    elseif sum(isnan(V))>0
        disp('NaN appeared');
        break
    else
        V00 = V;
        disp(['Iteration ',num2str(iter)])
        C0 = C;
        Eta0 = Eta;
        Beta0 = Beta;
    end
       
    iter = iter+1;
    
    toc
    
    disp(' ')
    disp(' ')
    
end

end
